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System identification or mathematical modeling is used in the aerospace community for development of simulation 
models for robust control law design. These models are often described as linear time-invariant processes. 
Nevertheless, it is well known that the underlying process is often nonlinear. The reason for using a linear approach 
has been due to the lack of a proper set of tools for the identification of nonlinear systems. Over the past several 
decades, the controls and biomedical communities have made great advances in developing tools for the identification 
of nonlinear systems. These approaches are robust and readily applicable to aerospace systems. In this paper, we 
show the application of one such nonlinear system identification technique, structure detection, for the analysis of F- 
15B Quiet Spike™ aeroservoelastic flight-test data. Structure detection is concerned with the selection of a subset of 
candidate terms that best describe the observed output. This is a necessary procedure to compute an efficient system 
description that may afford greater insight into the functionality of the system or a simpler controller design. 
Structure computation as a tool for black-box modeling may be of critical importance for the development of robust 
parsimonious models for the flight-test community. Moreover, this approach may lead to efficient strategies for rapid 
envelope expansion, which may save significant development time and costs. The objectives of this study are to 
demonstrate via analysis of F-15B Quiet Spike aeroservoelastic flight-test data for several flight conditions that 

1) linear models are inefficient for modeling aeroservoelastic data, 2) nonlinear identification provides a 
parsimonious model description while providing a high percent fit for cross-validated data, and 3) the model 
structure and parameters vary as the flight condition is altered. 


I. Introduction 

S YSTEM identification, or black-box modeling, is a critical step 
in aircraft development, analysis, and validation for flight 
worthiness. The development and testing of aircraft typically takes 
many years and requires considerable expenditure of limited 
resources. One reason for lengthy development time/costs is 
assuming that the underlying system is linear and invariant 
throughout the flight envelope. This assumption is related to having 
an inadequate knowledge of an appropriate model type or structure to 
use for parameter estimation. Selection of an insufficient model 
structure may lead to difficulties in parameter estimation, giving 
estimates with significant biases and/or large variances [1], This 
often complicates control synthesis or renders it infeasible. The 
power of using structure-detection techniques as a tool for model 
development (i.e., black-box modeling) is that it can provide a 
parsimonious system description that can describe complex 
aeroelastic behavior over a large operating range. Consequently, 
this provides models that can be more robust and therefore reduce 
development time. 

Moreover, when studying aeroelastic systems, it may not be 
practical to assume that the exact model structure is well known a 
priori. In aerospace systems analysis, one of the main objectives is 
not only to estimate system parameters but to gain insight into the 
structure of the underlying system. Therefore, structure computation 
is of significant relevance and importance to modeling and design of 
aircraft and aerospace vehicles. Structure computation may indicate 
deficiencies in an analytical model and could lead to improved 
modeling strategies and also provide a parsimonious black-box 
system description for control synthesis [1-3]. 
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For linear systems modeling, a commonly used approach for 
determining model structure is the minimum description length 
(MDL) proposed by Rissanen [4], MDL was specifically developed 
to overcome some of the inconsistencies of Akaike’s information 
criterion [5]; that is, its variance does not tend to zero for larger 
sample sizes N. 

Recently, the bootstrap method has been shown to be a useful tool 
for structure detection of nonlinear models [6-8]. The bootstrap is a 
numerical method for estimating parameter statistics that requires 
few assumptions [9], The conditions needed to apply bootstrap to 
least-squares estimation are quite mild: namely, that the errors be 
independent, identically distributed, and have zero mean. 

In this paper, we investigate whether 1 ) a linear or nonlinear model 
best represents the observed data and 2) the system is invariant during 
envelope expansion (varying Mach number). The data analyzed in 
this study are from the F-15B Quiet Spike™ flight-test program, 
which was a collaborative effort between Gulfstream Aerospace 
Corporation (Savannah, Georgia) and NASA Dryden Flight 
Research Center [10-12], The results show that 

1) Linear models are inefficient for modeling aeroservoelastic 
data. 

2) Nonlinear identification provides a parsimonious model de- 
scription while providing a high percent fit for cross-validated 
data. 

3) The model structure and parameters vary as the flight condition 
is altered. 


II. NARMAX Model Form 

The dynamic behavior of many nonlinear systems can be 
represented as a discrete polynomial that expands the present output 
value in terms of present and past values of the input signal and past 
values of the output signal [13-15]. A system modeled in this form is 
popularly known as a NARMAX (nonlinear autoregressive, moving- 
average exogenous) model and is linear in the parameters. 

Recently, Kukreja and Brenner [ 16 1 showed that NARMAX 
identification is well suited to describing aeroelastic phenomena. The 
NARMAX structure is a general parametric form for modeling 
nonlinear systems [L4], This structure describes both the stochastic 
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and deterministic components of nonlinear systems. Many nonlinear 
systems are a special case of the general NARMAX structure [17], In 
this paper, we focus on a special class of NARMAX models: 
nonlinear polynomial models. The polynomial NARMAX structure 
models the input-output relationship as a nonlinear difference 
equation of the form 

z(n) = f'[z(n — 1 z(n — n z ), u(n ), u(n — n u ), e(n 

-1 e(n - n e )\ + e(n) (1) 


[ j_,20,2J_] . In our treatment of the data, the ARM A and NARMA 
model classes are of practical significance. 

This special case of the general NARMAX model [Eq. (1)] can be 
written as 

z(n) = f'[z(n - 1), ... , z.(n - n z ), e(n - 1 e(n — «,)] 

+ e(n) (4) 

where we redefine the model order for this model set as 


where / denotes a nonlinear mapping, / is the order of the 
nonlinearity, u is the controlled or exogenous input, z is the measured 
output, and e is the uncontrolled input or innovation. This nonlinear 
mapping may include a variety of nonlinear terms, such as terms 
raised to an integer power, products of present and past inputs, past 
outputs, or cross terms. In general, the nonlinear mapping / can be 
described by a wide variety of nonlinear functions such as sigmoids 
or splines [ T7, j_8] .This system description encompasses many forms 
of linear and nonlinear difference equations that are linear in the 
parameters. 

Identifying a NARMAX model requires two things: 1) structure 
detection and 2) parameter estimation. Structure detection can be 
divided into model-order selection and selecting which parameters to 
include in the model. We consider model-order selection as part of 
structure detection because, theoretically, there are an infinite 
number of candidate terms that could be considered initially. 
Establishing the model order, then, limits the choice of terms to be 
considered. For the NARMAX model, the system order is defined to 
be an ordered tuple as 

0 = [n„ n z n e I] ( 2 ) 

where n u is the maximum lag on the input, n z the maximum lag on the 
output, n e the maximum lag on the error, and / is the maximum 
nonlinearity order. Note that for nonpolynomial NARMAX models, 
/ may be simply replaced by a nonlinear mapping of some specified 
class. In this paper, we assume that the system order is known. 

III. Structure Detection 

The structure-detection problem is that of selecting the subset of 
candidate terms that best describe the output. Therefore, the 
parametrization of a system is still further reduced by determining 
which of the components are required. A binomial coefficient is 
defined as 

/ n \ n(n — 1) • • • {n — k + 1) n\ 

\ k J k\ k\(n — k ) ! 


0 = [n z n e 1 ] (5) 

The maximum number of candidate terms in a model [Eq. (4)] with 
n z and n e dynamic terms and /th-order nonlinearity is 

( n z + n e + l\ ... 

P= I i J (6) 

Note that ARMA models can be estimated using the Yule-Walker 
equations or the instrumental variable estimator to avoid estimating 
the MA part [1_,20,2J_]. This is the approach taken in this paper. For a 
NARMA model, the NMA part must be modeled. For nonlinear 
systems, output additive noise can produce multiplicative terms 
between input, output, and itself. To compute unbiased parameter 
estimates, a noise model (i.e., NMA) needs to be estimated [22]. 

V. Structure-Detection Methods 

With the model types defined for the flight-test data available for 
analysis, we describe two approaches applicable to these model 
classes. The first is appropriate for autoregressive (AR) models and 
the second is appropriate for NARMA models. 

A. Minimum Description Length 

A commonly used technique in linear system identification to 
determine model structure is MDL [4], MDL was specifically 
developed to overcome some of the inconsistencies of Akaike’s 
information criterion [5]; that is, its variance does not tend to zero for 
larger sample sizes N. 

The number of parameters necessary to reproduce an observed 
sequence {zj . . . . , z n} of a time series depends on the model and 
parameters assumed to have generated the data [4], The MDL 
technique finds the model that minimizes the description length and 
thereby computes an estimate of model order [4]. 

Binary prefix codes are used to encode data strings. These data 
strings can be made up of symbols, parameters, numbers, etc. It is 
known that the average length of a code word is bounded by 
Shannon’s theorem [4]. Therefore, it is possible to write [4] 


where k\ is the factorial of k, k\ = 1 • 2 • 3 ■ . . . ■ (k — 1) • k. The 
maximum number of terms in a NARMAX model with n z , n u , and n e 
dynamic teims and /th-order nonlinearity is 


{ n u + n z + tz e 4- 1 + / \ 

V i ) 


(3) 


^j?0)L(x) > -^p(x)logp(x) (7) 


where L(x) is the length of the code word (i.e., length of parameter 
vector 9), and p(x) is the probability of x. It is also possible to write 

L(z\x, 9) = - log p(z\x, 9) (8) 


As a result, the number of candidate terms becomes very large for 
even moderately complex models, making structure detection 
difficult. We define the maximum number of terms p as the number 
of candidate terms to be initially considered for identification. 
Because of the excessive parameterization (the curse of 
dimensionality), the structure-detection problem often leads to 
computationally intractable combinatorial optimization problems. 

IV. Time Series 

The data considered in this paper are time series because the input 
signal u(n) is assumed to be zero or constant. Time-series analysis is 
used when inputs are not available to the experimenter or when it is 
unclear which signals are inputs and which are outputs [19]. Models 
arising from times-series data can have several unique forms 


where L(z \x, 9) is known as the log-likelihood function (to be 
maximized). Let 9 denote the value of the parameter that maximizes 
the likelihood and thus minimizes the parameter vector length (i.e., 
code-word length) L(z\x , 9). Because 9 can only be encoded up to a 
certain precision, the code-word length L(z\x, 9) becomes longer 
than the desired minimum L(z\x, 9), given noise considerations. Let 
the precision be S = 2~ q , where q is the number of bits used for 
encoding the parameter. It is possible to save on the code-word length 
if q is small. However, the result is a loss in precision. The optimal 
precision depends on the size of the observed data via — log <5 = 
0.5 log N and hence the total code-word length for k parameters is 
given by the MDL: 

MDL (k) = — log[maximized likelihood] + \k log N (9) 
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which, for an AR(/7 Z ) model, gives 

n 

MDL (n z ) = log[maximized likelihood] H — -log A? (10) 


B. Bootstrap 

Recently, the bootstrap has been shown to be a useful tool for 
structure detection of nonlinear models [8], The bootstrap is a 
numerical method for estimating parameter statistics that requires 
few assumptions [9], The conditions needed to apply bootstrap to 
least-squares estimation are quite mild: namely, that the errors be 
independent, identically distributed, and have zero mean. 

The bootstrap is a technique to randomly reassign observations 
that enables reestimates of parameters to be computed. This 
randomization and computation of parameters is done numerous 
times and treated as repeated experiments. In essence, the bootstrap 
simulates a Monte Carlo analysis. For structure computation, the 
bootstrap method is used to detect spurious parameters: those 
parameters for which the estimated values cannot be distinguished 
from zero. 

Application of an appropriate l 2 estimator to measured data gives 
the model response Z, residuals e, and parameter estimate 6. The 
bootstrap assumes that these residuals i = [e lf s 2 , . . . , %] arise from 
an unknown distribution V. By randomly resampling these residuals, 
with replacement, it is possible to generate a resampled version of the 
prediction errors e* = [s* , s|, . . . , ejy], for which the distribution 
estimates T>. The resampling procedure for each e* involves 
randomly selecting from e with an equal probability associated with 
each of the N elements. For example, a possible resampled version of 
the errors for N = 5 is e* = [e 4 , e t , e 4 , e 2 , e 3 ]. The star notation 
indicates that e* is not the calculated error e, but rather a resampled 
version of it. These resampled errors are added to the model response 
to generate a bootstrap replication of the original data: 

Z* = V z J + e* (11) 

A new bootstrapped parameter estimate 6* is obtained from this 
bootstrapped data Z * . This procedure is repeated B times to provide a 
set of parameter estimates from the B bootstrap replications: 

( 12 ) 

Parameter statistics can then be easily computed from 0* by fonning 
percentile intervals at a chosen level of significance a. 

Two methods related to the bootstrap are genetic algorithms and 
particle filters. They are both resampling procedures for which the 
goal is to estimate model parameters and distribution, mainly for the 
nonlinear and non-Gaussian case. 

Genetic algorithms are related to subset selection that may be 
loosely considered a type of Monte Carlo approach. However, as 
with subset selection, genetic algorithms fail to produce an optimal 
model set. Genetic algorithms are known to be computationally 
demanding, can lead to premature convergence on poor solutions, 
and have a tendency to converge toward local optima (and in some 
cases to arbitrary points) [23-26]. 

Particle filters implement a resampling method that tries to achieve 
the same goal as the genetic algorithm of estimating a sequence of 
hidden parameters based only on observed data. Particle filtering has 
been useful in a variety of areas but is difficult to implement for 
problems such as regressor selection, due to computational 
complexity leading to intractable problems for moderately complex 
full-model forms [27-301. 

Structure detection can provide useful process insights that can be 
used in subsequent development or refinement of physical models. 
Therefore, in the sequel, we investigate the applicability of MDL and 
the bootstrap to experimental aircraft data. Specifically, MDL and 
bootstrap methods are used as structure-detection tools to assess 
whether the 1) underlying data are best described by a linear time- 


invariant or nonlinear model and 2) model structure is invariant 
during envelope expansion. 

VI. Experimental F-15B Quiet Spike Data 

The MDL and bootstrap technique were applied to experimental 
flight-test data from the F-15B Quiet Spike project by Gulfstream 
Aerospace Corporation (Savannah, Georgia) and NASA Dryden 
Flight Research Center [10-12], The data analyzed for this study 
used the structural accelerometer response output of the Quiet Spike 
boom tip when fully extended. 

In this study, we only analyzed data from cases in which the Mach 
number was varied to be able to quantify whether model structure 
and parameters vary as the flight condition parameter is altered. If 
data were analyzed from cases in which more than one flight 
condition parameter was altered, one could not assess whether the 
dynamics and nonlinear effects vary due to Mach number or a 
combination of conditions. 

A. Procedures 

Flight data were gathered during subsonic, transonic, and 
supersonic flutter clearance of the F-15B Quiet Spike (see Fig. Q. 
This paper considers accelerometer data measured during pilot- 
induced pitch raps at Mach 0.85, 0.95, and 1.40 at 12,192 m (40, 
000 ft). The output was taken as the response of an accelerometer 
mounted near the Quiet Spike boom tip (see Fig. 2). Data were 
sampled at 400 Hz. For analysis, the recorded flight-test data were 
decimated by a factor of 2, resulting in a final sampling rate of 
200 Hz. 

For identification of a linear model, an arbitrarily large AR model 
of fiftieth order (n z = 50) was posed for estimation and the MDL 
technique was used to determine the optimal linear model. For 
identification of a nonlinear model, a model order with fourth-order 
output and error dynamics and second-order nonlinearity, 
O = [4, 4, 2], was used. A model with fourth-order dynamics was 
selected because it has been observed that aeroservoelastic structures 
are well defined by a fourth-order linear time-invariant system [31- 
33 ]. The nonlinearity order was chosen as second-order because 
empirical results showed that models of higher nonlinear order were 
not efficient to describe the data sets available for analysis. This gave 
a full-model description with 45 candidate terms. The nonlinear 
model was identified by applying the bootstrap approach. For the 
bootstrap method, B = 100 bootstrap replications were generated to 
assess the distribution of each parameter. For all three techniques, 
each parameter was tested for significance at the 95% confidence 
level. 

For both linear and nonlinear identification. Table 1 shows the 
number of data points available for each flight condition. The 
estimation data sets were from accelerometer response measure- 
ments of the primary sensor on the boom tip and the cross-validation 
data sets were of the backup sensor at the same location. 



Fig. 1 Flight-test article in extended configuration. 
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Primary 



B. Results 

The results of identifying the F-15B Quiet Spike data are 
presented. Figure 3 shows the output data sets used for this analysis. 
The data represent structural accelerometer response (primary sensor 
boom tip) used to compute the system structure. 

Equations (13-15) depict the model structure computed by the 
bootstrap method: 

Mach 0.85. altitude 40. 000 ft (12, 192 m), 
z(n) = 9 x z(n - 1) + 0 2 z(n - 2) + 0 3 z(n - 3) + 9 4 z 2 (n - 4) 

+ 9 5 e(n — 1) + d 6 s(n — 2) + 9 7 e(n — 3) 

+ 9 & z(n - 4 )e(n - 4) 4- 9 g s 2 (n - 4) + e(n) (13) 

Mach 0.95. altitude 40, 000 ft (12, 192 m), 
z(n) = 9 { z(n - 1) + i %z(n - 2) + ft 3 z(n - 3) 

+ 9 4 z(n - 3 )z(n - 4) + ft 5 z 2 (n - 4) + & 6 s(n - 1) 

+ t? 7 e(» — 2) + $ 8 e(n — 3) + z %z(n — 3 )e(n — 4) 

+ & 10 z(n - 4 )s(n - 3) + & n e{n - 3 )e(n - 4) 

+ ft l2 z{n - 4 )e(n - 4) + S l3 e 2 (n - 4) + s(n ) (14) 

Mach 1.40, altitude 40, 000 ft (12, 192 m), 
z(n) = Pi z(n - 1) + P 2 z(n - 2) + P 3 z(n - 3) 

+ P 4 z(n - 1 )z(n - 4) + P 5 s(n - 1) + P 6 e(n - 2) 

+ P 7 s(n - 3) + P%z{n — 1 )e(n - 4) + j) 9 z(n - 4 )e(n — 1) 

+ Pio s ( n ~ !)eO - 4) + e(n) (15) 

Equations ( 1 3—15) represent the computed model structure for flight 
conditions Mach 0.85, 0.95, and 1.40, respectively. The computed 
model structures are represented by a combination of linear and 
nonlinear lagged-output terms and contain 9, 13, and 10 terms for 
Mach 0.85, 0.95, and 1.40, respectively. Hence, the bootstrap 
technique successfully produced a parsimonious model description 
from the full set of 45 candidate terms. 

For AR (linear) model identification using MDL to compute 
structure, the estimated models were of order n, = 42, 44, and 46 for 
Mach 0.85, 0.95, and 1 .40, respectively. These models are not shown 
because they are simply a dynamic expansion of the output up to the 
order stated. However, for cross-validation data, we show the model 


Table 1 Data points available at each flight condition 


Altitude 12,192 m (40, 000 ft) 

Mach number 

0.85 

0.95 

1.40 

Estimation data 

572 

400 

504 

Validation data 

572 

400 

504 


fit of these linear models compared with the cross-validation fit 
obtained with the NARMA models (see Fig. 4). 

Figure 4 shows the predicted output for the cross-validation data 
sets for the identified structures; Fig. 4a displays Eq. ( 13 ) and 
AR(n z = 42), Fig. 4b displays Eq. (14) and AR(n z = 44), and 
Fig. 4c displays Eq. (15) and AR(n z = 46)). Each panel displays the 
full time history of the predicted output of the linear and nonlinear 
models superimposed on top of the measured output. For Mach 0.85 



a) 



b) 



c) 

Fig. 3 Estimation data: recorded structural accelerometer response to 
stick raps. 
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(Fig. 4a), the predicted outputs of the linear and nonlinear models 
account for over 91 and 95% of the measured output variance, 
respectively. For Mach 0.95 (Fig. 4b), the predicted output of the 
linear and nonlinear models account for over 92 and 97% of the 
measured output variance, respectively. For Mach 1 .40 (Fig. 4c), the 
predicted output of the linear and nonlinear models account for over 
91 and 96% of the measured output variance, respectively. 



a) 




c) 

Fig. 4 Cross-validation data: predicted linear and nonlinear model 
accelerometer response of z-tip longitudinal sensor superimposed on top 
of measured accelerometer output. 


The results demonstrate that although the AR models contain 
many more terms to explain the underlying process, they still offer a 
lower percent fit compared with the nonlinear model at the cost of 
model complexity (higher order), which often leads to more complex 
control synthesis. The nonlinear models contain only a few terms and 
were capable of explaining a larger percent of the output variance. 
For these data sets, the results show that linear models are inefficient 
for accurate modeling of aeroservoelastic data. These results show 
that a nonlinear identification approach offers a parsimonious system 
description while providing a high percent fit for cross-validated 
data. Moreover, the results illustrate the need to vary model structure 
for different flight conditions. 

VII. Discussion 

Experimental results demonstrate that structure computation as a 
tool for black-box modeling may be useful for the analysis of 
dynamic aircraft data. The bootstrap successfully reduced the 
number of regressors posed to aircraft aeroelastic data, yielding a 
parsimonious model structure for each data set. Additionally, these 
parsimonious structures were capable of predicting a large portion of 
the cross-validation data collected on a backup sensor at a similar 
location. However, for linear analysis, the MDL approach was not 
able to reduce the model order (structure) as well and yielded a more 
complex system description. Although these linear models have 
higher complexity (degrees of freedom), they provided a model- 
predictive capability that explained a smaller percent of the observed 
output variance. This find indicates that a linear model may not be 
appropriate to describe aeroservoelastic data. A higher percent fit 
offered by the parsimonious nonlinear models suggests that the 
identified structures and parameters explain the data well. Using 
percent fit alone as an indicator of model goodness could lead to 
incorrect interpretations of model validity. Nevertheless, in many 
cases, for nonlinear models, this may be the only indicator that is 
readily available. 

In this work, the results show that although the linear dynamics 
remained invariant for all flight conditions available for analysis, the 
nonlinear dynamics changed as the Mach number increased. For 
Mach 0.85, the model [Eq. (13)] displayed a mildly nonlinear 
process, which physically makes sense because the airflow is mainly 
subsonic. When the Mach number was increased to 0.95, the model 
[Eq. (14)1 demonstrated a richer nonlinear dynamic description, 
which is likely due to embedded shock formations in the transonic 
regime. For Mach 1.40, the model [Eq. (15)] displayed a mildly 
nonlinear process again, which physically makes sense, because in 
this regime, the shocks become fixed. It is difficult to make definitive 
comments on the underlying physics responsible for this behavior 
without extensive analysis of different flight conditions. The 
important points to note are that this study suggests that 

1) Nonlinear models are appropriate to describe the dynamics 
behavior of advanced aircraft. 

2) Models describing aircraft dynamics vary with flight condition. 

This suggests that nonlinear modeling may afford a robust and 

parsimonious system description compared with a larger operating 
regime and that models used for prediction (e.g., control) should not 
be invariant for all flight conditions. 

For this study, only a polynomial mapping with fourth-order 
output and error lag was used as a basis function to explain the 
nonlinear behavior of the F-15B Quiet Spike data. Clearly, different 
basis functions and a higher dynamic order (lag order) should be 
investigated to determine if another basis can produce accurate 
model predictions with reduced complexity. Moreover, further 
studies are necessary to evaluate whether the model structure is 
invariant under different operating conditions, such as altitude, and 
model parameterizations. 

This study illustrates the usefulness of structure detection as an 
approach to compute a parsimonious model of a highly complex 
nonlinear process, as demonstrated with experimental data of aircraft 
aeroelastic dynamics. Moreover, analysis of flight-test data can 
provide useful process insights that can be used in subsequent 
development or refinement of physical models. In particular, 
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morphological models are based on assumptions (e.g., these effects 
are important and those are negligible), which may be incorrect 
[34,35]. A structure computation approach to model identification 
may help uncover such surprises. 

VIII. Conclusions 

Results show that linear models are inefficient for modeling 
aeroservoelastic data and that nonlinear identification provides a 
parsimonious model description while providing a high percent fit 
for cross-validated data. Moreover, the results demonstrate that 
model structure and parameters vary as the flight condition varies. 
These results may have practical significance in the analysis of 
aircraft dynamics during envelope expansion and could lead to more 
efficient control strategies. In addition, this technique could allow 
greater insight into the functionality of various systems dynamics by 
providing a quantitative model that is easily interpretable. 
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